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. ABSTRACT 

\^ ' We present dynamical models of the nearby compact elliptical galaxy M32, using 

00 , high quality kinematical measurements, obtained with the integral-field spectrograph 

■ SAURDN mounted on the William Herschel Telescope on La Palma. We also include 

^-H ' STIS data obtained by Joseph et al. We find a best-fit black hole mass of M, = 

O ! (2.5 ± 0.5) X lOHf© and a stellar /-band mass-to-light ratio of (1.85 ± 0.15) Mq/Lq. 

For the first time, we are also able to constrain the inclination along which M32 

I is observed to 70° ± 5°. Assuming that M32 is indeed axisymmetric, the averaged 

r-j . observed fiattening of 0.73 then corresponds to an intrinsic fiattening of 0.68 ± 0.03. 

i' These tight constraints are mainly caused by the use of integral-field data. We 

JL I show this quantitatively by comparing with models that are constrained by multiple 

^ ■ slits only. We show the phase-space distribution and intrinsic velocity structure of the 

' best-fit model and investigate the effect of regularisation on the orbit distribution. 

I Key words: galaxies: elliptical and lenticular, cD - galaxies: kinematics and dynamics 

J> . - galaxies: structure, galaxies: individual, M32 - integral-field spectroscopy 



1 INTRODUCTION quality were built to match these data sets (Richstone, 

Bower & Dressier 1990; van der Marel et al. 1994b; Qian 



M32 is a high-surface brightness, compact E3 companion of 
the Andromeda galaxy. Ground-based kinematic measure- 



et al. 1995; Dehnen 1995). For most galaxies, the value of 

the best-fit black hole mass that is found depends on the 
ments or this galaxy (Tonry 1984, 1987) already showed a , i i ^ j_i i- j_ -i i- r i- 

, , It n, 1 assumptions that are made about the distribution function 

steep gradient m the central velocity pronle and a central /■rNTn\ r,i i m • . i i i i ■ i .i . 

(Dl ) oi the galaxy. Iwo-mtegral models, which assume that 

dispersion peak, suggesting the presence oi a central com- , , t-t^ i i i j_i ^ i • i ■ ^ i r i- 

i- I oo o f- ^Yie Dt depends only on the two classical integrals oi motion 
pact object, presumably a supermassive black hole. Since n j j_i j_- i j_ rj_i i 

, , , , (the energy E and the vertical component oi the angular mo- 

then, (spectroscopic) observations with increasing spatial , , \ n . i . i- . .i , i i i i 

' , , ,/T-, nx^i mentum Lz), generally tend to overpredict the central black 

resolution, both ground-based (Dressier & Richstone 1988; , , • j_i ^ • i re ■ ^ i- i i- 

° ^ hole mass, since they cannot provide sumcient radial motion 

van der Marel et al. 1994a; Bender, Kormendy & Dehnen , . , , ,„„„ ui ^ i r,r,nr, n ^ i 

, ^^„r^ ,^ , , , , , , (Magorrian et al. 1998; Gebhardt et al. 2000; Bower et al. 

1996) and with HST (Lauer et al. 1992; van der Marel, de 2mi) 

Zeeuw & Rix 1997; Joseph et al. 2001), strengthened the 

case for this black hole. rpjjg current state-of-the art models allow for the maxi- 

Simultaneously, dynamical models of ever improving ^.^^ ^^^^^^ anisotropy by assuming that the distribution 

function depends on three integrals of motion (van der Marel 
et al. 1998, hereafter vdM98; Cretton et al. 1999, hereafter 
* E-mail: verolme@strw.leidenuniv.nl C99; Gebhardt et al. 2000, 2001). The internal kinemati- 

t European Space Agency external fellow cal structure of the best-fitting three-integral model of M32 

t Hubble fellow closely resembles that of an f{E,Lz) model. This explains 
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why the central black hole mass that was found in M32 by 
the early models does not differ much from the current value 
of (3.4 ±0.7) X IQ^Mq (vdM98). 

This means that the mass of the central black hole and 
intrinsic kinematical structure of M32 are constrained rather 
well. The main uncertainty that remains is the inclination 
along which M32 is observed, or, equivalently, the intrinsic 
flattening: edge-on models to a data-set consisting of FOS 
kinematics and ground-based observations along four posi- 
tion angles fit equally well as models that are projected over 
55° (vdM98). When we combine this with the observed flat- 
tening of M32 (which is almost constant and equal to 0.73 
inside the central ten arcseconds, vdM98), we see that the 
intrinsic flattening of M32 is not very tightly constrained 
and can have any value between 0.55 and 0.73. 

Recent developments in instrument design offer a 
way out: integral- field spectrographs capture the full two- 
dimensional behaviour of objects and are therefore expected 
to better constrain parameters such as the inclination. This 
tighter constraint on the intrinsic parameters can be ex- 
plained from the behaviour of two-integral models, even 
though these generally provide less accurate fits to observa- 
tional data. The part of a two-integral distribution function 
that is even in the velocities is only determined once the 
meridional plane density distribution is known (Qian et al. 
1995). The full intrinsic velocity map is needed to constrain 
the odd part of the DF, which means we need to measure 
the line-of-sight velocity vioa- Whereas these two lowest or- 
der velocity moments are sufficient to determine the stellar 
DF of a galaxy, the dark matter distribution is only con- 
strained once the velocity anisotropy is known (Dejonghe 
1987; Gerhard 1991, 1993). This quantity can be measured 
through the higher order velocity moments. It is therefore 
necessary to determine the full two-dimensional kinematic 
behaviour of a galaxy. Capturing this information with only 
a few slit positions, as is generally attempted, is sometimes 
possible, but not always, which implies that the distribution 
function remains unconstrained. 

These effects are even more important for t/iree-integral 
distribution functions, which, in most cases, provide better 
fits to the observations and have more freedom. Therefore, 
when galaxy models are constrained by kinematic observa- 
tions along a modest number of slits, the intrinsic structure 
may remain largely unconstrained. In this paper, we revisit 
M32 using two-dimensional kinematical maps obtained with 
the panoramic integral-field spectrograph SAURON (Bacon et 
al. 2001, hereafter Paper I; de Zeeuw et al. 2002, hereafter 
Paper II) and high-resolution spectra obtained with STIS on 
board HST (Joseph et al. 2001). We show that the use of 
integral-field data places very tight constraints on the cen- 
tral black hole mass and mass-to-light ratio of M32, but 
also, for the first time, on the inclination, or equivalently, 
the intrinsic flattening. 

The paper is organized as follows: in Section ^, we give 
a brief summary of the data that are used to constrain the 
dynamical models, which are described in Section ^. The 
results are presented in Section ^ and we show in Section ^ 
that the use of integral-field data is crucial for most of the 
results that we obtained. The properties of the best-fitting 
model are described in Section ^ and we conclude with a 
discussion in Section |^. 

Throughout the paper we assume a distance of _D = 0.7 



Mpc (Welch et al. 1986). This choice does not influence 
our conclusions, but sets the scale of the models in phys- 
ical units. Specifically, lengths and masses scale as D, while 
mass-to- light ratios scale as . 



2 OBSERVATIONS 

We use two fully complementary data sets obtained with 
current state-of-the-art instruments: the integral-field spec- 
trograph SAURON mounted on WHT, and STIS on board 
HST. Integral- field spectrographs provide the full kinematic 
picture in one consistent data-set, while the high-resolution 
instrument STIS is ideal to probe the central kinematics of 
nearby objects. With this combined data-set, we expect to 
be able to determine the intrinsic properties of M32 with 
great accuracy. 

2.1 SAURON observations 

SAURON is mounted on the WHT on La Palma and is specifi- 
cally designed to measure the kinematics and line-strengths 
of nearby ellipticals, lenticulars and spiral bulges. Details on 
the instrument and the data reduction process can be found 
in Paper I. The instrument has two observing modes: in the 
low-resolution mode, the field of view is 30 x 41 arcsec at a 
pixel size of 0'.'94, while the high-resolution mode, used to 
zoom in on objects in good seeing conditions, has a pixel size 
of 0'.'28, corresponding to a field-of-view of 9 x 11 arcsec. 

The SAURON observations were made on October 15, 
1999 (Paper II). Relatively good seeing (FWHMsi 0'.'95) al- 
lowed us to obtain the measurements in the high-resolution 
mode. The kinematical maps of M32, derived from a single 
2700s exposure, are shown in the four upper panels of Fig. hi 
The velocity field reaches a peak value of « 65 kms~^; the 
misalignment between the rotation axis and the projected 
minor axis is less than two degrees, demonstrating that the 
measurements are consistent with axisymmetry. The disper- 
sion of M32 drops well below the instrumental resolution 
(cr ~ 90 kms~^) outside the inner few arcseconds. This re- 
sults in a fourth-order Gauss-Hermite moment hi (measur- 
ing the symmetric deviations from a Gaussian, van der Marel 
& Franx 1993) that is difficult to measure. As is apparent 
from the fourth panel of Fig. 0, it is almost constant over 
the field. 

A comparison of the SAURON measurements with 
ground-based long-slit data, obtained with the ISIS spec- 
trograph on WHT, shows that the velocities differ only by 
2.3 ± 2.2 kms~^, the dispersions by 6.0 ± 2.6 kms~^ and 
/i3 by —0.009 ± 0.016 (Paper II; the accuracy of /14 is diffi- 
cult to determine) . This indicates that there is no systematic 
offset in V and /13, while the disagreement between the dis- 
persion measurements is only significant at the 2cr-level. We 
do not correct for this and use all SAURON measurements to 
constrain the models. 

2.2 STIS observations 

The STIS measurements (Fig. ^) were obtained by Joseph 
et al. (2001), using the FCQ-method (Bender 1990). They 
clearly show the main signatures of a central black hole: a 
steep velocity gradient and a peak in the dispersion profile 
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Figure 1. Top. The kinematic maps of M32 (from left to right: the mean velocity, velocity dispersion and Gauss-Hermite moments 
/i3 and /i4), as observed with the integral-field spectrograph SAURDN. The field of view is 9 X 11 arcsec. Bottom. The best fit kinematic 
maps, obtained by adding the weighted contributions of the best fitting set of orbits in the SAURON apertures. 



near the center. By assuming that M32 has a two-integral 
distribution function and by solving the Jeans equations, 
Joseph et al. estimate that a central black hole with a mass 
of (2 — 4) X 10® Mq best reproduces these signatures, in full 
harmony with the earlier results by vdM98. STIS provides 
a continuous spatial sampling at « 0'.'051 pixels, avoiding 
positioning problems that may occur when using aperture 
spectrographs such as FOS (van der Marel, de Zeeuw & Rix 
1997). Comparisons between these two types of data are 
therefore only possible if such uncertainties are taken into 
account (Joseph et al. 2001). The STIS measurements are 
superior to the FOS data, since the spectral resolution and 
signal-to-noise ratio are much higher, an advantage when 
measuring the relatively low velocity and dispersion of M32. 
Because of this, we can use up to the fourth order Gauss- 
Hermite moment, while FOS was only able to measure the 
first two. 



3 METHOD 
3.1 Mass models 

We construct dynamical models, based on Schwarzschild's 
orbit superposition method (Schwarzschild 1979; 1982). We 



assume that the galaxy is axisymmetric and allow for the 
full degree of anisotropy by assuming that the distribution 
function depends on three integrals of motion. The imple- 
mentation that is used here was originally developed by Rix 
et al. (1997), vdM98 and C99, and was designed for use 
on one-dimensional surface-brightness profiles (but see Cret- 
ton & van den Bosch 1999). However, for M32, as for most 
nearby ellipticals, high-quality two-dimensional imaging is 
available, so that more sophisticated density parametrisa- 
tions are possible. An example of such a formalism is Multi- 
Gaussian Expansion (MGE, Monnet, Bacon & Emsellem 
1992; Emsellem, Monnet & Bacon 1994), in which the set 
of Gaussians that best fits an image of the galaxy is found. 
The advantages of this formalism are that both axisymmet- 
ric and triaxial galaxies can be reproduced in a straightfor- 
ward manner, the deprojection is simple and can be done 
analytically (under certain assumptions) once the viewing 
angles are known, and convolution with (MGE) point spread 
functions is fast and simple. We have therefore adapted the 
axisymmetric Schwarzschild software to use MGE; the new 
algorithm was tested while modeling the dynamics of the 
counter-rotating core galaxy IC1459 (Cappellari et al. 2002). 
These tests have shown that the modeling results do not de- 
pend on the choice of mass model: the best-fitting param- 
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Table 1. The parameters of the eleven Gaussians in the MGE-fit 
to the combined ground-based and HST/WFPC2/F814W image 
of M32. The second column shows the amplitude P', the third 
column shows the width a and the fourth the projected flattening 
of each Gaussian. 



Figure 2. The panels show the major axis velocity, dispersion, 
third and fourth order Gauss-Hermite moments of M32 obtained 
with the high-resolution spectrograph STIS on board HST (dia- 
monds, Joseph et al. 2001), together with the predictions at the 
STIS apertures of our best-fit model (solid line). 



eters and intrinsic (phase space) structure of models that 
start from an MGE parameterisation with constant axial 
ratio are identical to those of models that use a power-law 
density profile. 

An MGE parameterisation of M32 was obtained using 
the software of Cappellari (2002), by fitting to a combina- 
tion of a high-resolution HST/WFPC2/F814W image and 
a ground-based J-band image obtained by Peletier (1993; 
0'.'549 pixels"^ with a seeing of l'.'27) with the Isaac Newton 
Telescope on La Palma. Since we are constructing axisym- 
metric models, we do not allow for shifts in the position 
angles of the Gaussians and fix the centers of all individual 
Gaussians to the center of the galaxy. Each Gaussian there- 
fore has three free parameters: the amplitude P' (in units 
of Lq pc"'^), the width a' in pixels (or arcseconds) and the 
projected flattening q' . Eleven Gaussians were used in the 
fit, since, upon adding more Gaussians, the RMS error of 
the fit does not change by more than 1%. The parameters 
of these eleven Gaussians are listed in Table |l| this superpo- 
sition reproduces the M32 surface brightness with an RMS 
error of «2.2 per cent (Cappellari 2002, Figure 6). This is 
within the measurement errors. 



3.2 Axisymmetric three-integral Schwarzschild 
models 

We use the MGE-parameterisation to calculate axisymmet- 
ric three-integral models with the modified Schwarzschild 
method. Our models are different from those of previous 
authors: we use independent kinematic data sets and we 
parameterise the mass density differently. These differences 
provide a useful test: if we find best-fitting parameters that 



number 


P'(Lq pc-2) 


ct' (arcsec) 


9' 


1 


615412. 


0.0364 


0.740 


2 


517371. 


0.131 


0.809 


3 


324569. 


0.348 


0.775 


4 


175920. 


0.753 


0.732 


5 


54346 


1.51 


0.742 


6 


22279 


3.18 


0.708 


7 


11428 


6.05 


0.746 


8 


4037.5 


11.7 


0.683 


9 


2506.2 


19.2 


0.819 


10 


1142.4 


33.5 


0.820 


11 


226.4 


95.8 


0.827 



match the ones found by e.g. vdM98, it means that both 
approaches are reliable and the results are robust. 

The set of best-fitting Gaussians is deprojected by as- 
suming that the galaxy is axisymmetric and choosing a value 
of the inclination|^. The resulting intrinsic luminosity den- 
sity is multiplied by a constant mass-to-light ratio M/ L to 
obtain the intrinsic mass density. The stellar potential can 
then be calculated by applying Poisson's equation. A central 
supermassive black hole is mimicked by adding the poten- 
tial of a central point mass M, . In the combined stellar and 
black hole potential, a representative orbit library is found 
by sampling the two classical integrals of motion, the energy 
E and the vertical component of the angular momentum Lz , 
together with the third integral I^. 

Here, we only give a brief outline of the scheme that is 
used to sample these integrals; a more thorough description 
can be found in C99. The orbital energy is sampled through 
a logarithmic grid in the radius Rc of the circular orbit at 
energy E. This grid is chosen to include more than 99.99 per 
cent of the mass. While some of the orbits that make up the 
missing mass might contribute to the density in our radial 
range, this contribution will be smaller than 0.01 per cent 
and can therefore be safely ignored. 

At each energy, is sampled by a linear grid in the 
fractional angular momentum L^/ L^a^ £ ( — 1,1), with Lmax 
the maximum orbital angular momentum at energy E. This 
grid does not include L^/Lraa^ = {0,1} to avoid problems 
with the numerical integration. This is not a severe limita- 
tion, as the orbits that are missed in this manner carry only 
a small fraction of the mass in an axisymmetric model. The 
third integral is parameterised by the orbital starting point 
on the zero velocity curve. 

Every combination of the integrals (E, L^, I3) defines a 
separate stellar orbit. Since the third integral is generally 
not known analytically, the orbital properties can only be 
obtained by solving the equations of motion numerically. At 
every time step of this orbit integration, which is carried out 



3 The flattest Gaussian in the MGE superposition, in our case 
the eighth Gaussian, limits the range of possible inclinations to 
i > 47°. Since this range includes values of i that can be ruled 
out for other physical reasons (see below), this mathematical lim- 
itation is not a problem. 
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for a fixed number of characteristic orbital periods (usually 
chosen to be the period of the circular orbit at the given en- 
ergy), the orbital energy is randomized by drawing a number 
of arbitrary values from the energy bin around E, which are 
then assigned to the orbit. This 'dithering' of integral space 
results in smoother orbital superpositions, which generally 
provide better fits to the data (Zhao 1996; C99). 

The orbital superposition can be made even smoother 
by including the so-called two-integral components (TICs, 
C99; Verolme & de Zeeuw 2002) in the orbit library. These 
TICs are building blocks that obey only the two classical in- 
tegrals of motion E and Lz and implicitly include stochastic 
orbits at the given values of energy and angular momentum. 
Since they are smoother than most orbits, the orbital super- 
position also suffers less from discretisation when they are 
included in the library. Here, we calculated model fits to dif- 
ferent data-sets (both long-slit and integral-field data) with 
and without TICs in the orbit library, and found that the fit 
improves only marginally upon inclusion of the TICs (RMS 
error of the fit changes by less than 0.1%). On the other 
hand, calculating observables for the TICs at every aper- 
ture can be very time-consuming when many data points 
are included in the fit, as is the case when using integral- 
field data. To save computation time, we therefore decided 
not to include the TICs in the derivation of the best-fit pa- 
rameters. 

The observables of all orbits (and TICs) in the library 
are stored on grids that are adapted to the resolution of 
the observations. Using these, the orbital superposition that 
best fits the data is found by applying the non- negative least- 
squares (NNLS) routine of Lawson & Hanson (1974). The 
(unavoidable) discrepancy between model predictions and 
data, caused by noise in the data, discretisation effects and 
wrong choices of the model parameters, is expressed in a 
value for the x^y defined as 



E 



AD, 



(1) 



where Nd is the number of constraint points, Di is the obser- 
vational constraint at the i-th data point, D* is the model 
prediction at that point and ADi is the uncertainty that 
is associated with this value (usually the observational er- 
ror). By varying the inclination, central black hole mass and 
stellar mass-to-light ratio, we investigate which combination 
{M,, M/L,i) results in the overall smallest (the overall 
best fit to the data). The value of for ^ single model is of 
limited value, since the true number of degrees of freedom 
is generally not known, but the difference in x^ between a 
model and the overall minimum value, Ax^ = ~ Xmim is 
statistically meaningful (see Press et al. 1992). By using the 
central limit theorem, it can be shown that this interpre- 
tation is even valid when the errors ADi are not Gaussian 
distributed (van der Marel et al. 2000), which is important 
when applying the statistic to models of observational data. 
This means that we can assign the usual confidence levels 
to the Ax^ distribution and determine the probability that 
a given set of model parameters will occur. 

Introducing integral-field data in numerical models also 
means a considerable increase in the number of measure- 
ments that has to be reproduced: even when only one SAURON 
pointing is used, a fit to both the intrinsic and projected 



mass density and to all four Gauss-Hermite moments re- 
sults in approximately 8000 constraints. Standard mathe- 
matical practice in minimization problems is to use more 
free parameters than constraints, which means we would 
need ~ lO'* orbits to fit a SAURON field. On the other hand, 
the behaviour of models that are only constrained by obser- 
vational data is very similar to that of models with addi- 
tional regularisation constraints (R97; C99; vdM98; Fig. ^ 
of this paper), even though regularised models always have 
more constraints than orbits (§1.2). 

We conclude from this that it is not necessary to use 
more orbits than constraints to obtain accurate and unbi- 
ased fits to the data. In our case, we decided to use an orbit 
library with A''o = 20 x 7 x 7 orbits. Since the orbital Lz 
can be both positive or negative, depending on the chosen 
direction of rotation around the symmetry axis, this results 
in a library of 1960 orbits, implying that we have almost 
four constraints per orbit. 

When using integral-field data, the integration time of 
individual orbits may become extremely lengthy, as their 
properties have to be stored in many apertures. This is an- 
other reason not to use very large orbit libraries. Integrating 
all orbits in our library requires about eight hours of compu- 
tation time on a IGHz machine with 512Mb memory, while 
one NNLS fit takes approximately one hour on the same 
machine. 



4 RESULTS 

4.1 Best-fit intrinsic parameters 

We sample the model parameters (M,, M/L,i) on a grid of 
990 models, defined by fifteen values of the central black 
hole mass, spaced linearly between M, = and M, = 
4.5 X 10'' Mq, eleven values of M/L between 1.5 and 2.5, 
and inclinations i = 90°, 75°, 70°, 65°, 60° and 50°. We do 
not consider smaller values of the inclination, as this would 
make M32 intrinsically flatter than E7, i.e., flatter than any 
observed elliptical galaxy. Models with equal ratios between 
the mass-to-light ratio and the central black hole mass differ 
only by a scaling factor in the velocity. This means that, for 
every inclination, only fifteen orbit libraries have to be cal- 
culated, each with a different ratio M,/{M/L). Any point 
along a line of constant M, / [M / L) can then be obtained 
by scaling the velocity prior to the NNLS fit. As a result, 
one full x^ [M, , M/L) grid (at fixed inclination) takes about 
twelve days of computation time. 

For each of the 990 models that were obtained in this 
manner, we calculated the fit to the combined SAURON and 
STIS data. Contour plots of the values of x^iM., M/L,i), 
offset such that the overall minimum value is zero, are shown 
in Fig. ^ The contour levels correspond to the confidence 
levels of a Ax^ distribution with three degrees of freedom 
(with the thick contour in the figure indicating the 3cr confi- 
dence level), so that the volume in which the intrinsic param- 
eters of M32 are most likely located can be readily assessed. 
The shape of the contours is very regular, even though no 
smoothing has been applied. 

The best-fitting edge-on model has a black hole mass 
of (4.0 ± 0.5) X IO^Mq (3cr-error), which agrees very well 
with the vdM98's value of (3.4±0.7) x 10^ M©, even though 
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Figure 3. Ax'^ as a function of central black hole mass Al,, 
stellar mass-to-light ratio M/L and inclination i. The inclination 
varies from 90° (top panel) to 50° (bottom panel), and offset 
such that the overall minimum is zero. Contours are drawn at the 
formal confidence levels for a Ax^-distribution with three degrees 
of freedom, with the inner three contours corresponding to the la, 
la and Scr confidence levels (the thick contour indicates the 3cr 
level) ; subsequent contours correspond to a factor of two increase 
in Ax^. 



the observational data sets are independent. However, we 
see from Fig. ^ that models with an inclination in the range 
i = 70° ± 5° provide significantly better fits to the data than 
models with other inclinations, which are ruled out at high 
confidence. This narrow constraint on the inclination corre- 
sponds to an intrinsic flattening of approximately 0.68±0.03 
(given the observed flattening of 0.73). 

The overall best-fit model, located at i = 70°, has a 
black hole mass of (2.5 ±0.5) x 10^ and an 7-band mass- 
to-light ratio of (1.85 ± 0.15) Mq/Lq. The fit of this model 
to the SAURDN data is shown in the bottom panels of Fig. Q 
the fit to the STIS data is given in Fig. Q. 



4.2 Applying regularisation 

The matrix problem that is solved by the NNLS rou- 
tine is often ill-conditioned, resulting in orbital weights 
'y{E, Lz, I3) that vary rapidly as a function of the inte- 
grals of motion. This can be avoided by imposing regu- 
larisation constraints on the orbital weights (see e.g. Press 




100.0 



Figure 4. x as a function of the regularisation parameter, at 
fixed M, , M/L and inclination. The model with no regularisation 
is located at A = 00, but is placed at A = 500 to include the value 
in the plot. 



et al. 1992; Merritt 1993; Zhao 1996; Rix et al. 1997); 
another method, with similar effect, is to add maximum 
entropy constraints (Richstone & Tremaine 1988). Ap- 
plying regularisation forces the orbital weights towards a 
smooth function of E, and I3 by minimizing the n-th or- 
der derivatives d"^{E, L,,l3)/dE", d^-yiE, L,,h)/dL'i and 
d"'y{E, Lz, 73)79/3 . The degree of smoothing is determined 
by the order n and by the maximum fractional error A that 
the derivatives are allowed to have. For all regularised mod- 
els that are described here, we have minimized the second- 
order derivative (n = 2). 

The optimal value of the regularisation error A was de- 
termined by calculating model fits for a fixed combination 
of {M,, M/L,i) and varying A between 0.1 (high regulari- 
sation) and 00 (no regularisation). Fig. ^ shows that is 
large for small values of A, but quickly decreases to an al- 
most constant level for A > 4. Since the ideal model both 
fits the constraints and has a smooth distribution of orbital 
weights, we use A = 4 whenever regularisation is needed. A 
similar degree of smoothing was used earlier in applications 
that use the same regularisation technique (vdM98; Cretton 
& van den Bosch 1999) and in models that use the maximum 
entropy method (Gebhardt et al. 2000). 

Including smoothing constraints in the NNLS fit will, in 
general, change the value of x^ and therefore may affect the 
shape of the contours in plots such as Fig. |^. The top panel of 
Fig. ^shows a cross-section of Fig. ^at AI/L = 1.8; the of 
fits with the same parameters, but including regularisation, 
is shown in the middle panel. Since the shape of the contours 
does not change significantly, applying regularisation does 
not change the overall best-fit parameters, while it allows 
us to better study the phase-space structure and internal 
motions (see Fig. M below) . 

The two upper panels were calculated using four con- 
straints per orbit {Nd/No = 4, the value that is used 



3.2 



throughout the paper). Although we already argued in 
that our results are not sensitive to this parameter, we in- 
vestigated this more quantitatively by calculating (unregu- 
larised) models for which Nd/No = 0.6, which means there 
are almost two orbits per constraint point. The result of this 
is shown in the bottom panel, which indeed illustrates that 
our results do not depend on the number of orbits that is 
used. 
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Figure 5. Ax^ as a function of central black hole mass and incli- 
nation, for a fixed M/L of 1.8, for models without regularisation 
(top panel) and with regularisation (middle panel). The shape 
of the contours does not change upon adding regularisation. The 
bottom panel was calculated with N^/No = 0.6 and no regulari- 
sation, showing that our results do not depend on the number of 
orbits in the fit. 



5 THE EFFECT OF INTEGRAL-FIELD DATA 

At this point, although we have shown that the model pa- 
rameters are constrained much more tightly than in any pre- 
vious model, it is not entirely clear to what extent this is 
caused by the use of integral-field data. We also have to 
take into account that we have used a very dense sampling 
of M, and M/L, that we were able to calculate models at 
more inclinations than vdM98, and that the STIS data is of 
higher quality than the FOS observations. A more quantita- 
tive estimate of the influence of the integral-field data on 
the x'^ contours and on the best fitting parameters can be 
obtained by repeating the procedure of Section ^ for a data 
set that closely resembles the one used by vdM98. Since they 
fitted to a combination of FOS pointings with ground-based 
long-slit observations along the major, minor, 45° and 135° 
axes (which agree with SAURON), a similar data set can be 
obtained by extracting virtual 'slits' from the SAURON field 
along the same position angles and combining these with 
the STIS data. 

This data set was used to constrain models with the 
same parameter ranges as described in Section ^ Fig. ^ 
shows the of these fits. The best-fitting black hole mass 
is M, = (2.8 ± 1.2) X 10*^ Mq, the 7-band mass-to-light ra- 
tio is (1.75 ±0.25) Mq/Lq and the inclination can have any 
value below 75° . A comparison with the best-fit parameters 
of Section ^ indicates that the volume of intrinsic param- 
eters providing a reasonable fit to the data is significantly 
larger for the slit data set than when the entire SAURON field 
is used. We also notice that, although the use of integral- field 
data tightens the constraint on all model parameters (with 
a fractional change in the allowed range of Ad, of 2.4 when 




Figure 6. Similar to Fig. |3| for model fits to a data set consist- 
ing of the STIS data and virtual 'slits' simulated from the SAURON 
data. The range of best-fitting values in all three model param- 
eters is much larger than in Fig. ^, which indicates that using 
integral-field data tightens the constraint on all model parame- 
ters. 



changing from integral-field to slit data, and 1.7 for M/L), 
the most dramatic effect occurs for the inclination, which 
has an eight times larger allowed range when the galaxy is 
observed along a few slits only. 

This can be explained from the behaviour of two- 
integral models, even though these generally provide less 
accurate fits to observational data. As was mentioned in 
the Introduction, in order to constrain the even part of a 
two-integral distribution function, one needs to specify the 
meridional plane density distribution (Qian et al. 1995) . The 
odd part is determined once the full intrinsic velocity field is 
known, which can only be measured through the projected 
velocity field Vzi{x',y') (in which x' and y' are the coordi- 
nates in the plane of the sky). Furthermore, even if the total 
stellar DP is determined accurately, the dark matter distri- 
bution is still unconstrained. This can be solved by measur- 
ing the higher order velocity moments, which are related to 
the anisotropy of the velocity distribution. 

When the projected velocity moments are simple func- 
tions of x' and y' , measurements at a small number of posi- 
tions are sufficient to determine the behaviour at all interme- 
diate position angles. Indeed, some very specific two-integral 
models, such as the power-law galaxies (Evans & de Zeeuw 
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1994), have velocity fields that are quadratic in x' and y' , 
implying that observations along two perpendicular slits are 
sufficient to determine their full two-dimensional kinemati- 
cal behaviour. 

However, most galaxies require a distribution function 
that depends on three integrals of motion. The velocity fields 
of these three- integral models (and, in fact, of most generally 
apphcable two-integral models) are usually much more com- 
plicated functions of the projected coordinates (Dehnen & 
Gerhard 1993; de Bruijne, van der Marel & de Zeeuw 1996) . 
This means that these velocity fields cannot be character- 
ized by a (small) number of slit orientations. It is therefore 
not surprising that realistic models constrained by observa- 
tions along only a few slit orientations cannot provide much 
information about the intrinsic shape of the galaxy. It also 
implies that integral-field data places constraints on the in- 
trinsic shapes of galaxies, as well as on their mass distribu- 
tion. 

We see in both Fig. ^ and Fig. ^ that the best-fitting 
black hole mass of each separate panel increases with in- 
chnation. This behaviour can be understood at least par- 
tially from the following: the contribution of the intrinsic 
rotational velocity to the line-of-sight velocity decreases 
proportionally to sini. However, at smaller inclinations, the 
model must also be intrinsically fiatter, which means is 
larger when lowering i. For a two-integral model, the in- 
crease in is dominant, so that there is more net velocity 
at smaller values of the inclination. This also implies that 
less mass is needed to obtain the same velocity field. Unfor- 
tunately, since there is a trade-off between the mass-to-light 
ratio and the central black hole mass and because we are 
dealing with t/iree-integral models, it is not clear to what ex- 
tent these considerations can explain the effect we observe. 



6 PROPERTIES OF THE BEST-FIT MODEL 

The upper panels of Fig. ^ show the phase-space distribu- 
tion of the orbits of the non-regularised best-fit model as a 
function of /a and Lz, for constant E. The energy decreases 
from left to right, corresponding to an increasing distance 
from the galaxy center. The bottom panels are similar, for 
the model with a regularisation error of A = 4. 

We can study the intrinsic dynamical structure of 
the best-fitting model (M. = 2.5 x IO'^Mq, M/L=1.85 
Mq/Lqj, i = 70°) by adding the appropriate moments 
of the orbits with non-zero weight. This implies that we 
can study the degree of anisotropy of the model, and from 
this determine to what extent the galaxy resembles a two- 
integral model. Fig. ^ shows the intrinsic velocity moments 
{vl)^ (solid line), {vl)^ (dotted line) and {v^)^ (dashed 
line) of the best-fitting model, as a function of the meridional 
plane radius and the angle from the symmetry axis. The 
velocity moments are normalized by the total RMS motion 
Iv^) — {vl) -\- {vg) + {v^) . We see that the model resembles an 

f{E,Lz) model near the equatorial plane {{v^)^ = {vg)^) 
and becomes increasingly more anisotropic towards the sym- 
metry axis. This indicates that a two-integral model is not 
suitable to fit all constraints simultaneously, in agreement 
with the findings of vdM98 (their Fig. 10). Furthermore, as 
in their models, azimuthal motion dominates. 
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Figure 8. Three intrinsic velocity moments (f^)^ (solid line), 
(Dg)2 (dotted line) and {'v'^)'^ (dashed line) as a function of the 
meridional plane radius and the angle from the symmetry axis, 
for the best-fitting model without regularisation. The velocity 
moments are normalized by the total RMS motion, {'u^). 



For completeness, we have checked the internal kine- 
matical structure of the pseudo-slit model for the same pa- 
rameters as were used for Fig. ^ Since this combination of 
intrinsic parameters is inside the 3(t contour of Fig. ^, and 
therefore also provides a good fit to the data, the results re- 
semble those of Fig. ^ closely. The main difference between 
using slits and integral-field data is that more combinations 
of intrinsic parameters can be ruled out by integral-field 
data. A point (M, , M/ L, i) that is not inside the Scr contour 
of Fig. H, but cannot be ruled out by observations along a 
few slits, will therefore in general show a different intrinsic 
kinematical structure. 

The anisotropy of a model can also be expressed in 
terms of an anisotropy parameter. A possible definition of 
this parameter is 

Pb^I^ ^'^ ^"^r . (2) 

When the model is fully isotropic, Pb ~ 0, while /3b = 
1 for a model that consists entirely of radial orbits, and 
Pb = — oo describes a model with only circular orbits. Fig. ^ 
shows this parameter as a function of the meridional plane 
coordinates R and z, over the range that is constrained by 
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Figure 7. Best-fit model orbital weights as a function of the three integrals of motion E, and I3. Each panel shows the orbital weights 
as a function of I3 and Lz (in units of the angular momentum of the circular orbit; since orbits can move in both the prograde and 
retrograde direction, I/z/Lmax can be either positive or negative) at constant E. The energy decreases from left to right, corresponding 
to an increasing mean orbital radius (from O'.'Ol in the left panel to 10 arcsec in the right panel. Larger radii are not shown, since these 
arc not meaningfully constrained by the data). Top panels: the rapidly varying orbital weights for the unregularised model. Bottom 
panels: including a modest degree of regularisation (A = 4) smoothens the orbital weights, while still fitting the constraints. The color 
scheme is the same as in the upper panels. 
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Figure 9. A greyscale plot of the anisotropy parameter 
as a function of the meridional plane coordinates -R and z, to- 
gether with a typical isodensity contour (using i = 70°). The 
values were obtained by summing up the intrinsic moments of 
the set of orbits with non-zero weight in the best-fitting model. 
The darkest regions correspond to radial anisotropy (/3s > 0); 
the lightest region shows the points where the model is tangen- 
tially anisotropic (/3s < —1). In all other points, the model is 
only mildly anisotropic (— 1 < f^B < 0)- 



the SAURON data. A typical isodensity contour is overplotted 
(using i = 70°). We see that /3b is positive close to the 
galaxy center, so that the velocity distribution is radially 
anisotropic there. In all other parts of the meridional plane, 
the model is azimuthally anisotropic. This result does not 
change upon adding regularisation. 

Dynamical models of elliptical galaxies have revealed 
that a few other objects are also tangentially anisotropic 
(NGC4697, Dejonghe et al. 1996; NGC1700, Statler et al. 
1999), while the majority seems to be radially anisotropic 
(M87, Merritt & Oh 1997; NGC2324, R97; NGC6703, Ger- 
hard et al. 1998; M32, vdM98; NGC1600, Matthias & Ger- 
hard 1999; NGC3379, NGC3377, NGC4473, NGC5845, Geb- 
hardt et al. 2000; NGC1399, Saglia et al. 2000). However, 
this radial anisotropy is usually rather small and sometimes 
even accompanied by a transition to isotropy at small radii. 
Furthermore, not all models included the full velocity profile 
in the fit, which implies that the mass-anisotropy degeneracy 
may not have been broken. Finally, our models are the first 
to use fully two-dimensional kinematical observations, which 
makes comparison with earlier models even harder. There- 
fore, although our results seem to contradict these findings, 
we must also conclude that our knowledge of anisotropy is 
still limited. 
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r CONCLUSIONS 

We have presented dynamical models of the nearby com- 
pact elliptical M32, using data from the integral-field spec- 
trograph SAURON and from STIS on board HST. We have 
shown that our modeling software is able to deal with two- 
dimensional kinematical information and that the integral- 
field data tightens the constraint on all intrinsic model pa- 
rameters considerably. 

The axisymmetric three- integral model that best fits the 
data has a black hole mass of (2.5±0.5) x 10^ M© and a stel- 
lar /-band mass-to- light ratio of (1.85±0.15) Mq/Lq. These 
values confirm the best-fit parameters that were obtained by 
previous authors, although our modeling procedure is dif- 
ferent: we use a fully independent kinematical data-set and 
a different parameterisation for the mass density. Despite 
these differences in the methods, the best-fitting parame- 
ters that we find match the ones that were found by e.g. 
vdM98. This means that both approaches are reliable and 
that the results are robust. 

For the first time, we are able to determine the inclina- 
tion along which M32 is observed with great accuracy: the 
best-fit value is 70° ±5°. If M32 is indeed axisymmetric, the 
averaged observed fiattening of 0.73 then corresponds to an 
intrinsic flattening of 0.68 ± 0.03. We have shown in Section 
^ that this tight constraint is mainly caused by the use of 
integral-field data. This implies that integral-field data will 
provide us with more insight into the internal structure and 
kinematics of such objects. 

Although M32 is consistent with axisymmetry, it may 
be intrinsically triaxial, but seen along one of its princi- 
pal planes. Allowing for an intrinsically triaxial object also 
would enable us to study the effects of uncertainties in the 
deprojection on the best-fitting parameters more closely. An 
extension of our method to triaxiality is in progress (Verolme 
et al. 2002, in prep.), which will allow us to verify the as- 
sumptions made here. Furthermore, the tests that were car- 
ried out in this paper and the agreement with the results 
of other authors indicate that, given the assumptions, the 
results are robust. 

We have obtained observations with SAURON of a rep- 
resentative sample of nearby ellipticals, lenticulars and spi- 
ral bulges, for most of which high-resolution STIS data is, 
or will become, available. We are in the process of apply- 
ing the axisymmetric version of Schwarzschild's method on 
SAURON observations of a few of the sample galaxies that ap- 
pear consistent with axisymmetry (e.g. NGC821, NGC3377, 
NGC2974). The results of these dynamical models will pro- 
vide us with the intrinsic parameters of a considerable num- 
ber of objects, giving us unique insight into the formation 
and evolution of early-type galaxies. 
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